Neural signals regulating motor synchronization in the primate deep cerebellar nuclei

Movements synchronized with external rhythms are ubiquitous in our daily lives. Despite the involvement of the cerebellum, the underlying mechanism remains unclear. In monkeys performing synchronized saccades to periodically alternating visual stimuli, we found that neuronal activity in the cerebellar dentate nucleus correlated with the timing of the next saccade and the current temporal error. One-third of the neurons were active regardless of saccade direction and showed greater activity for synchronized than for reactive saccades. During the transition from reactive to predictive saccades in each trial, the activity of these neurons coincided with target onset, representing an internal model of rhythmic structure rather than a specific motor command. The behavioural changes induced by electrical stimulation were explained by activating different groups of neurons at various strengths, suggesting that the lateral cerebellum contains multiple functional modules for the acquisition of internal rhythms, predictive motor control, and error detection during synchronized movements.


Data analysis for Supplementary Figures and relevant text
Mutual information (related to Supplementary Fig. 3d and e).
In addition to the directional index (DI) and the prediction index (PI), the ability of each neuron to discriminate saccade direction (ipsilateral or contralateral) or task condition (synchronized or reactive) was also evaluated by calculating the mutual information (MI) 1,2 . For this analysis, neuronal activity was measured from the spike density function during 100 ms before (Unilateral and Bilateral neurons) or after (Postsaccade neurons) saccades. The MI was defined as, where P(c, r) was the joint probability of the task condition c and the neuronal response r (higher or lower than the median). For the comparison with the DI and PI, the sign of the MI was assigned so that the preference for ipsilateral saccades or synchronized task was positive and the preference for contralateral saccades or reactive task was negative.
As shown in Fig. 2b and c, both Unilateral and Bilateral neurons exhibited ramping activity before saccade execution, which might be temporally scaled for different SOAs. To quantify temporal scaling capability of each group of neurons, we compared (1) the timing of the peak activity, (2) the magnitude of the peak activity, and (3) the slope of the ramping profile measured from the spike density function aligned with saccades across the three SOA conditions. Each parameter was tested by repeated measures of ANOVA (P < 0.05). We also examined the slope of the time-scaled data normalized for the SOA (Supplementary Fig. 4c). The same analysis was also applied to the spike density functions (σ = 30 ms) obtained from individual trials ( Supplementary Fig. 4d) 3 . Fig. 9)

Decoding analysis (related to Supplementary
To evaluate the neuronal coding for saccade timing, population decoder was constructed for each group of neurons using the Neural Decoding toolbox for Matlab 4 . For comparison, we used the same set of the data as the partial correlation analysis in Figs. 3 and 4 (synchronized saccades in the preferred direction, 550 ms SOA). The decoder identified whether the ISI was shorter or longer than the median (timing of the next saccade, Supplementary Fig. 9a) or whether the saccade preceded or lagged the target onset (temporal error, Supplementary Fig.   9b). The classification accuracy was computed for every 100-ms interval (10 ms step) after saccades. Since we recorded only single neuron at a given moment, pseudo-population data for each neuron group were created by randomly resampling trials in the same condition. For each decoding run, the data was divided into five sets of trials, and the maximum correlation coefficient classifier was trained on the four sets and tested for the remaining. The mean population vector of normalized neuronal activity in the training set was used as a template for classification. The classifier calculated Pearson's correlation coefficient between the test data and the template, and the condition with the highest correlation was taken as the predicted result.
This procedure was repeated five times for cross-validation by swapping the training and test sets. Classification accuracy was defined as the average percentage of correct trials identified in the test set during cross-validation. These procedures were repeated 100 times with different pseudo-population data, and the mean values are reported in Supplementary Fig. 9. To evaluate whether the decoder performance was better than chance, the permutation test (100 repeats) was performed by shuffling the relationship between neuronal activity and saccade timing. The result of the permutation test for each time point (10 ms) is shown by the solid bar below the traces (P < 0.05). We also constructed the decoder for tempo ( Supplementary Fig. 9c, based on the data of synchronized saccades in the preferred direction with different SOAs), task condition ( Supplementary Fig. 9d, synchronized and reactive saccades in the preferred direction with different SOAs), and saccade direction ( Supplementary Fig. 9e, synchronized saccades in both directions with different SOAs) using the same classification method (maximum correlation coefficient classifier).

Comparison of neuronal activity related to sensory versus motor events.
Correlation analysis (related to Supplementary Fig. 10c) The neuron shown in Fig. 5b exhibited a strong transient activity that peaked just before synchronized saccades (> 5th saccades, bottom trace with orange shading), but the activity bottomed out around the time of the second saccade in the sequence (red trace with orange shading), indicating that the neuronal activity did not always coincide with saccades. In addition to the sensorimotor index (SMI) that compared the magnitude of neuronal activity for sensory versus motor alignments, we also quantified the similarity of the time courses of neuronal activity between the second and the later saccades by calculating the correlation of the spike density profiles (± 275 ms aligned either with the target onset or saccades). The correlation of the saccade-aligned data for Bilateral neurons was significantly smaller than that for Unilateral neurons, while the correlation of the target-aligned data showed an opposite trend. These results further support that the activity of Bilateral neurons was related to sensory prediction rather than saccade execution. Fig. 11a)

Time warping analysis (related to Supplementary
To further clarify whether the neuronal activity was better aligned to sensory prediction or motor execution, we performed the time warping analysis modified from Perez et al. (2013) 5 .
We used only the data during the transition from reactive to synchronized saccades (a 700-ms period starting from 200 ms before the second target onset). To quantify the trial-by-trial variability in neuronal activity associated with target onset (Ls), we calculated the likelihood of spike occurrence on each trial using the normalized spike density function obtained from the remaining trials as a probability distribution. Next, the data were transformed to be aligned with saccades. Saccade timing in each trial was time-warped to the mean of all trials, and spike trains before and after the saccade was accordingly stretched or shrank in time without changing the total length of the data (700 ms). Then, we computed the likelihood of spike occurrence associated with motor execution (Lm). The relative magnitude of the trial-by-trial variability for the sensory and motor alignments was evaluated by calculating the time-warping index (TWI) defined as, If the neuronal activity was better aligned with the target onset, the TWI for early reactive saccades would have a positive value. For comparison, we also calculated the TWI for synchronized saccades (a 700-ms period starting from 350 ms before the sixth and subsequent target onset). These data are shown in Supplementary 11a in comparison with the sensorimotor index (SMI) calculated for the second saccades. Supplementary Fig. 11b)

ROC analysis (related to
The activity of saccade-related neurons initially lagged behind the target onset but later preceded it as the trial progressed (Fig. 5a). To examine the changes in the time course of neuronal activity for the second versus later targets (> 5th), we quantified the difference in phase of the spike density functions aligned with the target onset using the receiver operating characteristics (ROC) analysis. The ROC curve was constructed by comparing the normalized spike density profiles during the first and the later cycles for different levels of threshold that spanned the entire time range (1,100 ms). The area under the ROC curve (AUC) was greater than 0.5 if the neuronal activity in the first cycle lagged behind that in later cycles, and less than 0.5 otherwise. The mutual information is signed so that the preference for ipsilateral saccades is negative. These two indices significantly correlated (r = 0.90, P = 1.8 × 10 −34 ). As with the DI, the mutual information did not differ between Bilateral and Postsaccade neurons (Tukey test, P = 0.56), while mutual information for Unilateral neurons differed from that for the others (P < 0.01). (e) Comparison between the mutual information for task condition and the PI. The mutual information is signed so that the preference for the reactive saccade task is negative. These two indices were significantly correlated (r = 0.37, P = 0.0003). As in the PI, only the mutual information for Bilateral neurons was significantly different from zero (one sample t-test, two-sided, t 30 = 3.36, P = 0.002). Counts of neurons with significant changes either in the slope of ramping activity, the timing or magnitude of the peak activity. Most Bilateral (29/33) and Unilateral (32/33) neurons changed either of the three parameters (ANOVA, P < 0.05). For almost all neurons, these changes were scaled orderly except for a few neurons exhibiting a selective change only for the 550-ms SOA (one Unilateral neuron for peak timing, one Unilateral and two Bilateral neurons for peak magnitude).   Fig. 3c-e). (b) Population activity and the partial correlation coefficient for the ISI calculated for the data aligned with the next saccades. Since the peak of the preparatory activity aligned with the next saccades was similar between trials with early and late saccades, the partial correlation coefficients around the time of saccade (± 100 ms, gray shading) were not different from zero (one sample t-test, two-sided, P > 0.05). There were positive partial correlations at earlier time periods for both Unilateral and Bilateral neurons, reflecting the difference in ramping activity between trials with short versus long ISIs. Note that the sign of the partial correlation is reversed due to the difference in alignment because the faster increase in preparatory activity in trials with shorter ISI (shown in a) corresponds the shorter lead time of the ramping activity (in b).
(c) The time course of the partial correlation shown in a are extended backward in time, showing that there was no significant correlation between the neuronal activity and the two later saccade timing. In all panels, data are presented as mean ± s.e. The convention of the graph is the same as that in Supplementary Fig. 6d. (a) Classification accuracy for the timing of the next saccade by constructing population decoder for each neuron group, using the same data set for the partial correlation analysis in Fig. 3e (synchronized saccade task, preferred direction, 550 ms SOA). The coloured bars below the traces indicate when the decoder classified trials significantly above chance (permutation test, P < 0.05, maximum correlation coefficient classifier). Decoder constructed from Bilateral and Unilateral neurons successfully discriminated the timing of saccades, consistent with the results of the partial correlation analysis shown in Fig. 3e. Decoder for Postsaccade neurons also showed good results, probably due to the inverse correlation of successive saccade timing during synchronization (which was controlled in the partial correlation analysis). (b) Performance of the decoder for temporal error using the same data set as the partial correlation analysis in Fig. 4e (synchronized saccade task, preferred direction, 550 ms SOA). (c) Performance of the decoder for tempo (SOA) around the time of synchronized saccades in the preferred direction. (d) Performance of the decoder for task condition constructed from the data for both the synchronized and reactive saccade tasks. (e) Performance of the decoder for saccade direction constructed from the data for synchronized saccades. Saccade direction could be decoded reliably from the population activity of all types of neurons, consistent with the ANOVA results shown in Supplementary Fig. 2b. were computed for the 2nd through the 5th saccades in the preferred direction (upper panels). The same data were separated into four groups according to temporal error (saccade latency), and the SMI for each group was also computed (lower panels). When the target order was early and the temporal error was large (i.e., reactive saccades), the activity of Bilateral neurons (red lines) tended to be more aligned with the target than saccade, resulting in larger SMI values. When the target order was late and the temporal error was small (synchronized saccades), the SMI value became close to zero, indicating that the magnitude of neuronal activity was comparable between different alignments of the data. (b) Distributions of the SMI measured for saccades in the non-preferred direction for Bilateral and Postsaccade neurons. (c) Left: Distribution of the rank correlation coefficients between the spike density profiles (± 275 ms) for early reactive (2nd saccade) and later synchronized saccades (> 5th) aligned with saccade. Right: Distribution of the rank correlation coefficients calculated for the same data but aligned with the target onset. For Unilateral neurons, the rank correlation coefficients were large for the data aligned with saccades but were small for the data aligned with the target onset, indicating that they exhibited activity associated with saccades irrespective of the order of saccade in the sequence. By contrast, for Bilateral neurons, the correlation for the data aligned with saccades was the smallest among the three types of neurons but that for the data aligned with the target onset was the largest. These coefficient values were closely correlated with the SMI measured for the second saccades (motor aligned, r = −0.68, P = 1.5 × 10 −13 ; sensory aligned, r = 0.74, P = 6.6 × 10 −17 ). , during (stim), and after (stim+1) electrical stimulation. The box-whisker plot indicates the median, quartiles, and range of the data. Circle indicates the mean value. The residuals for ISIstim were significantly larger than those for ISIstim−1 (Tukey test, P = 0.009, n = 40), while there was no significant difference between the residuals for ISIstim−1 and ISIstim+1 (P = 0.09). These results suggest that electrical stimulation directly changed the timing of the next saccade, but the apparent stimulation effect on the two later saccade was mostly due to the behavioural adjustment during synchronization. Negative values indicate that the greater the neuronal activity, the more the next saccade is promoted (or the shorter the ISI). Since most Unilateral neurons showed a preference for ipsilateral saccades (Fig. 2a), correlations of Unilateral neurons with different directional preferences were calculated for the preferred (non-preferred) directions and shown for ipsilateral (contralateral) saccades. (c) Estimated impacts of neuronal activity on subsequence saccade timing at the time of electrical stimulation. Pearson's correlation coefficients were computed at different stimulation times (100 ms duration) following synchronized saccades. For the analysis in Fig. 6e, we assumed that the effects of electrical stimulation could be accounted for by the combination of signals in different type of neurons (Methods). (d) Distribution of the coefficient of determination (r 2 ) for fitting the data of electrical microstimulation. The black curve represents the null distribution generated by the permutation method (1,000 repeats). (e) Cumulative distributions of stimulation sites along the dorsoventral axis for the three groups. The colours match the groupings in Fig. 6c-e. (f) Effect of electrical stimulation separated for dorsal and ventral stimulation sites. The box-whisker plot indicates the median, quartiles, and range of the effect size. Stimulation of dorsal sites slowed, but stimulation of ventral sites facilitated, the next saccade. Negative values indicate that the greater the neuronal activity, the more the following saccade is promoted. Right: Correlation coefficients calculated for the data aligned with saccades. (b) Normalized activity of Bilateral neurons in one-third of the trials with early and late saccades in the reactive saccade task (mean ± s.e.). Note that these neurons exhibited ramping activity for contralateral reactive saccades only, while they showed preparatory activity for synchronized saccades in both directions ( Fig. 2a and Supplementary Fig. 3). (c) Distribution of the directional index (DI) of Bilateral neurons computed for synchronized (red) and reactive (black) saccades. The DI in the reactive saccade task was significantly smaller than zero (one sample t-test, two-sided, t 30 = 3.56, P = 0.001), and was significantly different from that in the synchronized saccade task (paired t-test, two-sided, t 30 = 3.51, P = 0.001).